Spatial distribution functions of random packed granular spheres obtained by direct 

particle imaging 
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We measure the two-point density correlations and Voronoi cell distributions of cyclically sheared 
granular spheres obtained with a fluorescence technique and compare them with random packing 
of frictionless spheres. We find that the radial distribution function g(r) is captured by the Percus- 
Yevick equation for initial volume fraction cf> = 0.59. However, small but systematic deviations are 
observed because of the splitting of the second peak as cj> is increased towards random close packing. 
The distribution of the Voronoi free volumes deviates from postulated F distributions, and the 
orientational order metric Qe shows disorder compared to numerical results reported for frictionless 
spheres. Overall, these measures show significant similarity of random packing of granular and 
frictionless spheres, but some systematic differences as well. 

PACS numbers: 45.70.Qj, 05.65.+b 
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The packing of spheres is one of the enduring prob- 
lems in physics, and a basis to understand the structure 
and strength of granular matter. Assuming dominance of 
steric interactions, dense packing of steel spheres was first 
used to understand structure of simple liquids with the 
radial distribution function g(r) and the orientation order 
metric Qq [ll. However, experimental measurements at 
boundaries [2j and computer simulations in the bulk Q 
have since shown that inter-particle friction can affect 
granular packing. Friction between particles changes the 
fundamental stability condition at contact from the fric- 
tionless case, causing a packing to be protocol dependent 
and the system to be out-of-cquilibrium. 

The difficulty of accurately measuring significant num- 
ber of particle positions in the bulk away from the influ- 
ence of boundaries has also stymied progress. Recent ex- 
perimental studies [4, 5] have examined packing of gran- 
ular spheres and find that the associated free volume dis- 
tributions are described by a T distribution with two fit- 
ting parameters which were then given a thermodynamic 
interpretation [5]. These results are puzzling in light of 
earlier analytical work in one-dimension and simulations 
in two and three dimensions that show a T distribution 
with 3 fitting parameters is needed to describe a broad 
range of volume fraction for elastic particles [fj. 

Here, we discuss new experiments with spherical gran- 
ular particles which enable us to directly determine sta- 
tistical measures to understand the effect of friction, test 
the effect of shear, and perform a rigorous comparison 
with frictionless hard sphere packing. Using a fluores- 
cence technique 0, 0, [1] , we obtain the packing of glass 
spheres before and after application of cyclic shear, and 
compare with random packing of frictionless spheres. We 
find that the overall shape of g(r) for volume fraction 
4> ~ 0.6 is captured by the Percus-Yevick equation 
which assumes random packing of spheres without angu- 
lar correlations. But, systematic deviations are observed 
because of the splitting of the second peak as <j> is in- 
creased toward random close packing, Qg shows partial 
hexagonal order, and the distribution of the Voronoi free 




FIG. 1. (a) Schematic diagram of the cyclic shear cell used in 
the experiments, (b) An image of the initial packing observed 
in the central vertical slice of the shear cell after particles 
are filled inside the cell, (c) Ordering grows near the top 
boundary which is free to move after 600 shear cycles but 
the packing in the bulk appears random, (d) The probability 
distribution function of the diameter of the glass beads, (e) 
The volume fraction cj> as a function of shear cycle number 
measured in the bulk (red/grey) evolves more slowly than 
in the entire cell (black) because of the ordering near the 
boundaries. 



volumes shows enhanced probabilities at higher values 
compared with T distributions postulated @ for random 
packing of spheres. 

The experiments to measure structure of granular 
packing are performed using a shear cell shown schemat- 
ically in Fig. HJa) . The parallelepiped shaped cell con- 
sists of a rigid front, back, and bottom transparent glass 
boundary, and side boundaries that can be tilted through 
a prescribed angle 6 to cyclically shear and perturb the 
packing. Glass spheres with density p g = 2.5 x 10 3 
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kg m~ 3 and average diameter d = 1.034 mm are gently 
added inside a shear cell filled with an interstitial liquid 
with the same refractive index as the glass spheres, den- 
sity pi = 1.0 x 10 3 kg m~ 3 , and viscosity v = 2.2 x 10~ 2 
Pas. Then a flat plate is placed on the top, which is 
constrained to move only in the vertical and horizontal 
direction and not allowed to rotate using a rigid set of lin- 
ear guides. The initial volume fraction of the glass beads 
is measured using the mass of the particles added and the 
volume of the cell occupied and found to be <f> — 0.59. A 
normal stress a z — 0.4 Pa is applied on the top boundary 
which is about five times the net gravitational stress due 
to the weight of the grains alone inside the cell, and is 
found to eliminate gradients due to gravity in the system. 

A dye is added to the liquid and a thin slice of the 
cell is illuminated with a laser and a cylindrical lens [13] • 
The resulting fluorescent light causes the particles to ap- 
pear dark against a bright background, and is imaged 
from an orthogonal direction with a resolution of 20 pix- 
els to a particle diameter using a 1000 x 1000 pixel 10- 
bit camera. A stack of images is recorded by linearly 
translating the plane of illumination. We then examine 
a 40c? x 5d x 17. bd central region as in indicated by box 
in Fig. [TJb,c) to avoid any effect of the boundaries, and 
locate the absolute position of the spheres to within the 
slight polydispersity of the particles (see Fig. QJd)). 

We impose quasi-static shear strain by varying 9 be- 
tween ±7r/36 radians with a mean angular speed u) = 
8.0 x 10 _3 rad s _1 , with a wait time of 50 s while the stack 
of images is acquired every time the system returns to its 
original position, 9q = Orad. The lubrication forces [ll[ 
due to liquid draining at contacts can be estimated to 
be 10 -5 lower than the confining forces, and the particle 
Reynolds number is ~ 10 _1 . Therefore the particles can 
be assumed to be in contact during the entire experiment 
and the interstitial liquid does not have any impact on 
the observed structure. 

The mean packing volume fraction of the spheres in- 
side the entire cell is first obtained by measuring the po- 
sition of the top plate as a function of the shear cycles. 
The mean 4> is observed to increase logarithmically from 
the initial value by 5% (Fig. dfe) consistent with previous 
reports with a similar setup [I4]. However, it is note- 
worthy that this is the total (j> inside the cell and can 
be different than <fi in the bulk because of influence of 
boundaries (Io| . Examining the images corresponding to 
the initial state of the packing, before applying the shear 
deformation, N = 1 and after shear cycle N = 600 shown 
in Fig. HJb,c), we indeed note greater ordering near the 
top where the boundary shears the spheres and moves to 
accommodate changes in the total <fi. The boundary be- 
tween ordered and disordered region appears sharp and 
moves downward as N is increased further, similar to de- 
velopment of crystalline order inside a Couette shear cell 
upon extended shearing 17]. Therefore, we focus on the 
first 600 shear cycles where the particles inside the bulk 
in the observation window appear uniformly random and 
obtain (f> from the ratio of the particle volume and the 




FIG. 2. The radial distribution function g(r) plotted as func- 
tion of distance r normalized by the mean diameter d for 
initial packing cf> = 0.59 and final packing (f> = 0.605 obtained 
after shearing, (black), is compared with the theoretical calcu- 
lation (red/grey) obtained by using the Percus-Yevick equa- 
tion, and the measured polydispersity of the beads. The 
4> = 0.605 case is offset for clarity. Inset: The calculated 
pair distribution functions Qij for particles coarsened to three 
sizes (4> = 0.59). The thick red/grey curve represents the av- 
erage of the six distinct contributions g%j{r) and is used for 
comparison with the experimental g(r). 

average Voronoi volume in the bulk. The Voronoi vol- 
ume is defined by points in the volume closest to that 
particle, and is calculated using algorithms written by 
Rycroft [l3| and measured particle positions. As shown 
in Fig. [He), the evolution of <fi in the bulk is observed to 
be slower compared with 4> measured for the entire cell 
and is used in all subsequent discussion. 

To analyze the structure from the measured particle 
positions, we first discuss the radial distribution function 
g(r) which represents the probability that the center of 
a particle is found at a distance r from another particle. 
g(r) obtained from the experimental data is shown in 
Fig. [2] For the initial volume fraction <p = 0.59, g(r) 
shows a tall peak at r ~ d, and a broad peak at r ~ 
2d, but for the higher <f> obtained after cyclic shear, the 
second peak splits and a weak secondary peak occurs at 
r = \/3d, corresponding to the next nearest neighbor in 
a face-centered-cubic (FCC) lattice. 

Now, the Percus-Yevick equation [j| [l4[ can be used to 
analytically find g(r) for randomly distributed particles 
for a given particle size, and volume fraction. Because of 
the slight polydispersity of the particles used in our ex- 
periments, we in fact have to compute the Percus-Yevick 
pair distribution function for polydisperse spherical pack- 
ings gij (r) , where the indices i, j represent the probability 
of finding a particle with diameter dj at a distance r from 
a particle with diameter d,. 
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FIG. 3. (a) The probability distribution of Q§ measured for 
each particle using Voronoi neighbors are observed to de- 
scribed by Gaussian fits, (b) The mean Qe,iocai measured 
as a function of (j>. The curve is a guide to the eye and shows 
that Qe increase somewhat over the narrow <j> investigated. 



here n is the number of particle sizes in the system. 
Coarse graining the observed size distribution - shown in 
Fig. Old) - into three sizes d\ = 0.98 mm, e?2 = 1.02 mm, 
rf 3 = 1.07 mm (using a greater n does not change the 
results significantly), we calculate the corresponding six 
distinct gij{r) terms, which are plotted in the Inset to 
Fig. [2] for (f> = 0.59, and thus the computed gpy(r) for 
the polydisperse packing according to the Percus-Yevick 
approximation, which is plotted in Fig. [5J We observe 
that the amplitude and width of the primary peak and 
the broad features of the secondary peaks are in good 
agreement with the Percus-Yevick approximation. This 
comparison is especially noteworthy because there are 
no fitting parameters. At higher <f>, the primary peak 
and the overall form is still captured by the ffpy(r), but 
details such as the splitting of the second peak which 
can indicate hexagonal ordering is not captured because 
Percus-Yevick approximation assumes random angular 
orientation. 

To investigate orientational order in the packing, we 
use the orientation order metric 
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where, I — 6 to examine hexagonal order, Yj m are the 
spherical harmonics, 0(r) and $(r) are the polar and az- 
imuthal angle, respectively, and r is the position vector 
from a particle to its neighbor [l5[. We define particle 
neighbors as those which share a Voronoi cell surface [l3[ . 
This removes any ambigity as is introduced when consid- 
ering only neighbors at contact due to roundoff errors 
in finding a particle center. In order to compare with 
packing of elastic particles, we first compute Qe, global 
by averaging Y/ m (0(r), $(r)) over all the bonds of the 
packing, and find Qe, g iobal — 0.27 ± 0.02. If particles 
neighbors are uncorrelated, then Qq is small because it 
goes as square root of the total number of bonds [iH ]. 
On the other hand, Qq for a FCC crystal with 12 neigh- 
bors is 0.5745. But even a slight perturbation due to 
roundoff errors introduces 2 extra neighbors and the cor- 
responding Q 6 is on average 0.454 [13] . Therefore, the 
observed distribution shows ordering but the degree of 
order appears lower compared with simulations of fric- 
tionless hard spheres [l8[, where Qq as high as 0.4 was 
reported for a frictionless hard spheres at comparable </>. 
In those studies which were performed with considerably 
smaller system size, particle inelasticity was observed to 
lower Qq but not as significantly as in our experiments. 

To examine the local orientational order more closely, 
we plot the observed probability distribution of Q§ for 
each particle in Fig. [3Ja), and the mean of the distribu- 
tion (Qe^ocai) in Fig. |5fb). (Qe^ oca i} is more sensitive 
than Qe^giobai to small crystalline regions within a pack- 
ing and allows us to avoid the possibility of destructive 
interference between different crystalline regions [li| . No 
significant enhancement of distribution is found at the 
values corresponding to FCC crystal, and the observed 
Qq distribution can be described rather well in fact by 
Gaussian fits. From these observations we conclude that 
while there is some local hexagonal order which increases 
slightly over the <p investigated, no significant crystal- 
lites occur in this dense regime approaching random close 
packing. 

A complementary method to examine the packing at 
the particle scale is using the free volume Vt associ- 
ated with each particle given by subtracting the min- 
imum Voronoi volume corresponding to close packing, 
v c — d 3 /^/2 from the Voronoi volume. This statistical 
quantity has gained prominence because it may be used 
to define a new measure of entropy based on disorder 
in packing @, |2jj, and may be amenable to thermody- 
namic interpretation It has been postulated based 
on analytical work in 1-dimensional systems, that Vf dis- 
tribution of random packing of spheres can be described 
by a r distribution 
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with three fitting parameters to, 5 and a that control dif- 
ferent parts of the distribution and were determined by 
numerical simulations with frictionless hard spheres [6(. 
In Fig. 21 we plot v / normalized by the mean free volume 
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FIG. 4. The probability distribution function of the free vol- 
ume associated with a sphere Vf normalized by the mean free 
volume (vf) plotted for various <f>. The smooth black curve 
is obtained after averaging over all the experimental data, 
r function corresponding to elastic frictionless spheres, the 
smooth red/grey curve, is shown for comparison, and is ob- 
served to deviate systematically at higher Vf. Allowing the 
fitting parameters to float improves the fit, but systematic 
differences persist for Vf > (vf) (blue dashed curve). 

(vf) at that </> along with the mean distribution obtained 
after averaging over all the measured <f>. The errors asso- 
ciated with the slightly polydispersity and errors in find- 
ing particle centers is of order of symbol size. Further, we 
plot Eq. [3] using m = 15.5, <5 = 1.3 reported in Ref. @], 
and allowing a to float to obtain best fit. Clearly, sys- 
tematic deviations are observed from the frictionless case, 
complementing the results for Qq. 



In order to check if a T-distribution can capture the 
experimentally observed free volume distributions, we 
tested both the three fitting parameter distribution, and 
the two fitting parameter distribution, which corresponds 
to setting 5 = 1 in Eq. [31 The best fit obtained with 
m = 12.3, a = 24.5, and S = 0.73 is also shown in Fig. [4] 
Even in this case we obtain enhanced probabilities for Vf 
greater than the mean. Therefore, our distribution differ 
from the experimental distributions used to give a sim- 
ple thermodynamic interpretation of granular packing 0] • 
While it is possible that such deviations arise because of 
the differences in preparation protocol, we note no signif- 
icant differences in the distributions obtained before and 
after application of cyclic shear in our experiments. 

In conclusion, we measured packing of granular spheres 
and compared experimentally obtained two-point density 
correlations and free volume distributions before and af- 
ter application of shear. We find that the radial distri- 
bution function is captured overall by the Percus-Yevick 
equation, which is important because it is fundamental 
to calculating the strength, heat conduction, and electro- 
magnetic wave scattering properties of a material. How- 
ever, angular correlation can be observed using the ori- 
entational order metric. In comparing with numerical 
simulations reported for frictionless sphere at compara- 
ble volume fractions, we find systematic differences in 
packing as measured by lower angular correlations, and 
deviation of free volume distributions from T distribu- 
tions postulated for frictionless hard spheres. 
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